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Abstract 

Very few works exist to date on development of a consistent energy-based coupling of atom- 
istic and continuum models of materials in more than one dimension. The difficulty in con- 
structing such a coupling consists in defining a coupled energy whose minimizers are free from 
uncontrollable errors on the atomistic/continuum interface. In this paper a consistent coupling 
in three dimensions is proposed. The main challenge in three dimensions was to identify and 
efficiently treat the modified Cauchy-Born continuum model which can be coupled to the exact 
atomistic model. The convergence and stability of the method is confirmed with numerical tests. 

Keywords: Consistent energy-based atomistic/continuum coupling, quasicontinuum method, multiscale 
method, three dimensions 

AMS subject classification: 65N30, 70C20, 74G15, 74G65 

1 Introduction 

Modeling defects in crystalline solids requires using atomistic models. On the one hand, defects 
create long-range elastic fields, accurate resolution of which requires a huge number of atomistic 
degrees of freedom, often unhandlable even on modern computers. On the other hand, the elastic 
deformation far away from a defect is well described by a continuum model. This is a rationale 
for atomistic/continuum (A/C) coupled methods — the methods that use full atomistic resolution 
around a defect, coupled to a coarse-grained continuum model far from the defect [TU t ITT ] IT8]. 

Consider a problem of finding an equilibrium of a certain atomistic crystal with a localized 
defect, i.e., finding a critical point of a potential energy of such a crystal. An A/C coupling 
approach to this problem would be to consider the exact energies of the atomistic deformation 
near the defect and a Cauchy-Born continuum energy of a Pi finite element discretization of the 
deformation field far from the defect. The efficiency of an A/C coupling rests on the fact that 
the complexity of computing the energy and the effective forces associated with an element is 
independent of the size of the element (which would not be true if the full atomistic model was 
used everywhere). 

*The work was performed during the author's stay at the Chair of Computational Mathematics and Numerical 
Analysis (ANMC) at the Swiss Federal Institute of Technology (EPFL) whose support is acknowledged. 
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The two main variants of an A/C coupling are the energy-based and the force-based coupling, the 
first one defines an A/C coupled energy that depends on the atomistic and continuum deformation, 
while the second one mixes the atomistic and the continuum forces (i.e., derivatives of the energy 
of the atomistic model and the continuum model); see PT8] and references therein for more 

details. The force-based coupling can indeed approximate effectively the exact atomistic equilibrium 
equations, however its stability properties are not well understood at present [TJ [21 [TT] and indeed 
there seem to exist examples of a force-based coupling of stable atomistic and continuum equations 
being unstable [T5] . 

When using energy-based methods, one faces a different kind of challenge: it turns out to 
be difficult to design a coupling which is at least first order consistent (a first order consistency is 
equivalent to a first-order convergence provided stability Q Despite the recent progress in designing 
a consistent energy-based A/C coupling [31 [TJ [HI [161 El) no satisfactory solution exists in the 
general case to date. 

One of the recent developments is the work [16], where the author proposed a consistent A/C 
coupling for two-body interaction in two dimensions. The key instrument in constructing a consis- 
tent coupling of [16] was the two-dimensional bond density lemma, which asserts that the effective 
number of atomistic bonds in a certain direction r G Z 2 lying on any triangle with vertices restricted 
to the lattice Z 2 equals to the area of the triangle, regardless of the direction r. This lemma al- 
lows to define the A/C coupling method in terms of the energy of individual bonds and show that 
continuum approximations of bond energies sum up to a (discretized) Cauchy-Born energy, up to 
some correction near the interface. 

The purpose of this paper is to extend the method of [16| to the three-dimensional case. Unfor- 
tunately, the three-dimensional analog of the bond-density lemma is not true: the number of bonds 
lying in a tetrahedron depends on the bond direction and in general is not equal to the volume of 
the tetrahedron. This makes the extension to 3D not trivial. 

The construction of the method in 3D is similar to the lower dimensional construction: we 
first define the continuum energy of bonds consistent with the exact energy of the bonds, and 
then show that the sum of continuum energies of bonds can be computed efficiently. The resulting 
three-dimensional continuum model turns out to be different from the Cauchy-Born model (this 
is a consequence of lack of the 3D bond density lemma). Nonetheless, numerical tests conducted 
confirm a similar behavior as the in 2D |13} [T6] . 

The paper is organized as follow. We formulate the proposed A/C coupling in Section [2] and 
define the effective number of bonds within a tetrahedron, BondVol(T, r); efficient computation 
of this quantity is central to the overall efficiency of the proposed method. Section [3] is entirely 
dedicated to an efficient algorithm of computing BondVol(T, r) and a Matlab code of this algorithm 
is given in Appendix [B] In Section [4] we present numerical tests of accuracy and stability, and in 
Section [5] we give concluding remarks. 

2 Consistent Atomistic/ Continuum Coupling 

For generality, we present the atomistic/continuum coupling in R rf , but will mainly focus on case 
d = 3. The cases d = 2 (considered in [16]) and d = 1 (considered in [HUB]) will be particular cases 

1 in the engineering-oriented literature, a lack of consistency is formulated in terms of fictitious forces called "ghost 
forces" 
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Figure 2.1: Geometry of an A/C interface. Black atoms belong to the (discrete) atomistic region 
f2 a , white atoms belong to the continuum region O c . The small atoms, belonging to Od, are involved 
in posing the Dirichlet-type boundary conditions. The discrete domains are, respectively, £ a , £ c , 
and £d- 



of the below discussion. 



In Section 2.1 we define the continuum and discrete regions, the triangulation of the continuum 



region, and the corresponding functional spaces. In Section 2.2 we present an atomistic interaction 



in terms of bond energies. Finally, in Section 2.3 we formulate the proposed A/C coupling. 
2.1 A/C Coupling Geometry 

Consider an atomistic material occupying a bounded domain Q C M d in its undeformed (reference) 
state. Assume a splitting of into three (open) regions, fi a where the exact atomistic model will be 
used, Q c where a continuum approximation will be used, and Qj) containing atoms whose positions 



will be fixed (when posing Dirichlet-type boundary conditions); see Fig. 2.1 for an illustration. The 
atom positions in the undeformed state comprise the lattice £ = nnZ d (where • denotes a closure 
of a set). We also denote £d = £ n S)d, £ c = (£ H Q c ) \ £d, and £ a = £ \ (£d U £ c ). Normally, 
the number of atoms in £ a is much less than the total number of atoms. 

It should be stressed that although we assume perfect crystalline lattices without defects in the 
reference configuration, the computational method can be presented if defects are allowed in the 
atomistic region. 

We define the space of all deformations, U, as all discrete functions £ — > M d and the space of 
admissible displacements Uq := {u G U : u\c D = 0}. Additional assumptions on Cd will be made 
later to avoid unnecessary complications due to boundary effects. 

We assume that f2 c is a polytope (i.e., polyhedron for d = 3) and Th is its triangulation with 
simplices T € Th- The spaces of A/C deformations and admissible displacements are defined as 

U h = {u h : O c U £ a — > M rf : u h is continuous on tt c and is affine on each T E 7~h}, 

and U% := {u h G U h : u h \c D = 0}. 
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2.2 Bond Formulation of the Atomistic Model 

We assume that the atomistic interaction is given by a set of neighbors 1Z C Z d \{0} and a two-body 
potential 4> : R d — > R. Define an interval (x, x + r) between two points x, x + r G M 2 as a set 

(x, x + r) := {x + Ar : AG (0, 1)} 

and call it a bond; r will be called a direction of a bond. Introduce the finite difference associated 
with a bond b = (x, x + r) or a bond direction r: 

Ay := A- y(x) := y(x + r) - y(x), 
and let energy of the atomistic model be 

E(y) :=£>(Al/)- (2.1) 
beB 

The variational equilibrium condition for y G IA is then 

(«E(l/),«) = VueW 

where (•, •) is the scalar product in and 5E(yp) £ U is the Gateaux derivative of E defined as 
(5E(y),v}:=£E( y + av)\ a=0 . 

We assume that dist(SO, Q \ Qt>) > max rS 7^ \r\ so that the following discrete version of the 
divergence theorem holds: 

D r v(x) = Vf G W , Vr G K. 

x<E£n(C-r) 

It is then easy to verify (see |16l Eq. (2.6)]) that a uniform deformation y^{x) := Fx is an equilibrium; 
i.e., (SE(vf),v) = 0. 

In the next subsection we propose an A/C couping E h (y), which is an approximation of E(y), 
such that 

6E h (y f ) = 0; (2.2) 

in other words, E h (y) is patch-test consistent. This method will be a generalization of the one- 
dimensional method [HI [16] and the two-dimensional method [16] (more precisely, the version of the 
method labeled as ECC in [T6]). 

2.3 The Proposed A/C Coupling 

Define the set of continuum bonds B c := {b G B : d C Sl c }, and the continuum directional 
derivative (associated with the bond direction r) 

V r y(x) := lim (f(x + er) - f(x)) . 

e-»0 

Further, define averaging over a bond b = (x, x + r): 

/•x+r r pi 

I f( x ) db = I f(x) db := / /(x + Ar)dA. 
A Jo 
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The proposed (patch-test-) consistent A/C coupling then reads 

E h (y) := <KA?/) + E / ^ V n>y) db > (2-3) 

where ,8 a := B\B C and r;, denotes the direction of b € Z3. The consistency (condition ( |2.2[ )) for this 
coupling follows from the fact that the variation of the continuum energy of a bond b S B c is equal 
to the variation of the exact energy on a uniform deformation yp: 

(6(J b <f>(V rb yF))db,v) = J b ^r b )V rb vdb = cj>{¥r b )D n v = (6<j>(D rb y F ),v) 

and its proof is completely analogous to the proof of the 2D version of this statement |16l Prop. 
3.2]. 

A straightforward calculation (see Proposition [T] below) allows us to convert bond integrals in 
(2.3) into a sum over elements (i.e., effectively volume integrals): 

E h (y) = <KA2/) + E E ^><KV r y| T ), (2-4) 

feeB a TET h r&Tl 

where the effective volumes of T are defined as 

/x+r 
X T db (2.5) 



xez d 



and the characteristic function x» is defined in Section 2.3.1 



It should be noted that the second term in (2.4) differs from the standard Cauchy-Born energy 
(for it to be equal to the Cauchy-Born energy, we must have fif^ = |T|). This in particular means 
that the existing results on stability [H [6] and consistency [T2] of the Cauchy-Born model does 
not apply to the coupling proposed in the present paper. Therefore we confirm convergence and 
stability of the proposed coupling via numerical tests in Section |4j 

2.3.1 Characteristic Function 

For a polytope u> C M. d (e.g., polyhedron in 3D) define a characteristic function 

Xu (x :=hm ' P \ M , 2.6 
p^o \Bp(x)\ 

where B p {x) is the ball centered at x with the radius p. We note that (i) the limit w.r.t. p — > in 
the definition of X w (aO exists, and (ii) including/excluding the boundary of a polytope oj (or any 
part of it) does not change the point values of x w { x )- The characteristic function is defined so that 
if UJ = UJi U UJ2 and uj\ fl u)2 = then Xuj = X u + X^ 2 pointwise. In particular, we have 

*n>) = S x t^ for a11 x G r2 - ( 2 - 7 ) 

TeT h 
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The characteristic function of a 3D polyhedron oj can be visualized as 



1 if x £ interior of uj 

2 if x G face of to 

if x € edge of w with angle a 
^- if x is a vertex of a; with spherical angle /3 
otherwise. 



(2. 



Note that the values of X ( x ) at the vertices of uj will not be important for the formulation of the 
method. 

With this definition of the characteristic function, we can prove the following proposition. 



Proposition 1. The energy (2.3) is equivalently written as (2.4) 



Proof. Indeed, the second sum in (2.3) can be transformed as 



W<KV r ,y)db= £ Yl f +T m-y)dh 

b£B c Jb veil x& d J * 



ren x& d 

(i,i+r)cS! c 



EE/ x nc mv)dh 

ren x& d Jx 

(i,x+r)C!i c 

E E f +r E x T m r y)db 

ren x& d J x TeTh 

(i,i+r)c!i c 

£X>(V r y| T ) /*T db > 



TeT h ren 



xei 

(x,x+r)e£l c 



where we used (2.7) and the fact that x a = 1 on any bond which lies inside J7 C . 



□ 



2.3.2 Complexity of Computing E h 



The method (2.3) with precomputing Qt,t directly using (2.5) may already yield a significant 
reduction in number of operations. Indeed, one must spend 0(#(£> c )) operations (here #(•) denotes 
the number of elements in a set) on precomputing Qtt only once for a given A/C geometry, and 
it would then take 0(#(i3 a )) + 0(#(7ft)#(7?.)) operations for computing the forces or assembling 



the stiffness matrix corresponding to (2.3) (recall that #(£> a ) ^ #(£>c))- 

Furthermore, in the ID case and in the 2D case with triangulation nodes coinciding with the 
lattice sites, the bond-density lemma yields that ^T,r = \T\ if T is far enough from the a/c interface 
and thus E c (y) reduces to the standard Cauchy-Born energy up to an interface correction. This 
removes the need in the precomputation step and yields an algorithm with an optimal complexity 
0(#(B a )) + 0(#(T h )#(TZ)). 

Unfortunately, as also shown in [16j . in general £It,t 7^ \T\ in 3D. (Numerical calculations of £It,v 
with randomly generated T and r show that £It,t 7^ \T\ for most choices of T and r.) Nevertheless, 
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Reduction to: 


# (parameters) 


BondVol(T, r) 


15 


gcd(ri, r 2 ,r 3 ) = 1 


15 


r = e 3 


12 


truncated prism 


9 


triangle 


6 


trapezoid 


4 


right triangle 


2 


S a ,b 


2 



Table 1: List of reductions from the original problem of computing BondVol(T, r) to computing 
S at b, with the number of parameters left after the reduction. 



as will be shown in the present paper, one can come up with an algorithm of computing Oj> 
efficiently. 

To this end denote, for any polytope u), 

rx+r 

BondVol(w,r) := ^ + x^db (2.9) 
so that £It,t can be expressed through BondVol(T, r) and the interface correction: 



Q T , r = BondVol(T, r) - ^ j- x T db. (2.J.U) 

(i,i+r)£B a X 

Here the quantity BondVol(w, r) is an effective number (volume) of all the bonds with direction r 
that intersect a polytope u. In the next section we will show that BondVol(T, r) can be computed 
with 0(log(diam(T)) +log |r|) arithmetic operations in 3D. This implies an overall precomputation 
time of at most o((#T h )(#K) (log(diam(ft)) + log(diam(7£))) . 

We expect that the precomputation time will not overly dominate the main computation time 
in most of applications. Indeed, the factor C(log(diam(f2)) + log(diam(7£))) is essentially the 
maximal number of iterations of the Euclid's algorithm and is between 15 and 50 for a typical 
atomistic system with diam(7£) ~ 5 and 10 2 S diam($7) S 10 9 . Indeed, in the numerical tests 



conducted in this paper, computing the volumes (the first term in (2.10)) was several times faster 



than doing the interface correction (the second term in (2.10)). 



3 Computing BondVol(T, r) 

This section is devoted entirely to an algorithm of fast computation of BondVol(T, r) defined by 



(2.9). A reader who is not interested in details or justification of the algorithm can skip this section 
or refer to Appendix [B] for a Matlab code. 

The algorithm is based on a series of steps, presented in Sections |3.1f 3.3 that reduce the 



original problem of computing BondVol(T, r) with 15 scalar parameters (12 to define T and 3 to 
define r) to an integer sum S a ^ (cf. ( |3. 13 ) ) with only two parameters a, b G Z. The principle steps 
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of the algorithm are summarized in Table [T] The overall complexity of the algorithm is discussed 
in Section l3T4l 

In this section by A, B, C, D, X, Y, Z we will denote the points in R 3 , and by r,s,x,£ we will 
denote vectors in M 3 . The points may be identified with their radius vectors. For two points 
X,Y 6R 3 , we denote XY = Y — X. For points and vectors, the subscripts 1,2, and 3 will denote 
their coordinates and we will use the notation X = (X\,X2, X%). The standard basis of M 3 will be 
denoted as e%, e%, e^. The points on the xy-plane in M 3 (i.e., the plane {X G M 3 : X3 = 0}) will be 
identified with the respective points in M 2 . 

3.1 Reductions 

We start with a tetrahedron T and a vector r = (n, r^i rs) G Z 3 , r ^ 0. 

3.1.1 Reduction to the case gcd(ri, ri-, r 3 ) = 1 
In this paragraph we show that 

BondVol(T,r) = BondVol(T, r/ gcd(n, r 2 , r 3 )), (3.1) 

where gcd(ri, r2, r 3 ) is the greatest common divisor of r\, r%, r 3 G Z. 

Indeed, let r = ns with n = gcd(ri, r2, r 3 ) G N and s G Z 3 . Fix x G Z 3 and consider a collection 
of points 

ATj. = {x + is : i G Z}. (3.2) 

First, compute the contribution of a collection of the respective collection of bonds {(£, £ + s) : £ G 
^} to BondVol(T,s): 

zZf *T db = E/ X T dh = z2 X T (x + (i + X)s)d\ = x T (x + Xs)dX. (3.3) 

Second, compute the contributions of 

+r) : CG^}= {(x + is,x + fs + r) : i G Z} 

= {(x + j s + ir, x + j s + ir + r) : i G Z, j = 0, 1, . . . , n — 1} 



to BondVol(T,r): 



2 f x T db = zZ/Zf X T dh 



n-l 



" poo 

^2 / X T (^ + js + Ar)dA 

1=0^°° ' (3.4) 



"X! / X T (^ + ^)d/x 



OO 



X T (x + ns)d/i, 
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where we did the change of variable /i = j + An. 



From calculations (3.3) and (3.4) it is easy to see that BondVol(T, s) = BondVol(T, r). Indeed, 



summing the contributions of different X x yields: 



BondVol(2» = E f +S *T dh ( 3 - 5 ) 
= E E *r db = BondVol(T, r), 

x x texjt 



which proves (3.1). 



3.1.2 Reduction to the case r = e 3 

We now assume gcd(ri, r 2 , r 3 ) = 1. In this subsection we first find a suitable linear transformation 
M such that Mr = e 3 , and apply it to both T and r. We then extend the definitions of BondVol(w, r) 
and Xuj t° allow measuring angles of edges and vertices in the untransformed space. 
Construction of M. Due to Bezout lemma, there exist ci 2 , d\2, c 3 , d 3 E Z such that 

ric 12 + r 2 di2 = gcd(ri,r 2 ), gcd (n, r 2 )c 3 + r 3 d 3 = gcd(gcd(ri, r 2 ), r 3 ) = 1. 

Take matrix M E Z 3x3 as the following product of two matrices: 



M 





Cl2 C?12 
gcd(ri,r 2 ) gcd(ri,r 2 ) 





and compute Mr: 



C12 ^12 0\ /rA /gcd(ri,r 2 ) 







t-2 n 
gcd(n,r2) gcd(n,r2) 

1/ \r3/ \ r 3 











)-l 





= e 3 









Mr = 

It is also important to notice that both M and M _1 have integer coefficients, the latter is due 
to det M = 1. Hence 

MZ 3 = Z 3 . (3.6) 

Extension of Xw and BondVol(w, r). We need to apply the transformation M to both T and r. 
To that end, we extend the definition of Xw by allowing for measuring angles of edges and vertices 
of u after applying M: 

X^!(x):=X M - 1 JM- 1 x), 
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so that Xu}(x) = Xm^(Mx). In case if oj is a polyhedron, x™} can be evaluated as 



l 

2 

Q 

2tt 
ii 





if x G interior of oj 
if x G face of ui 

if x G e, e is an edge of w, a is the angle of M _1 e in M _1 cj 
if x is a vertex of w, (3 is the spherical angle of M _1 x in M _1 u; 
otherwise. 



Likewise, extend 

BondVol(M, u,r) := BondVol(M _1 o;, M _1 r) 



(3.7) 



so that BondVol(w, r) = BondVol(M, Mo;, Mr). Notice that in the last step of (3.7) we used (3.6). 
The quantity BondVol(M, to, r) differs from BondVol(u;, r) only if r is parallel to some edge of uj 
and the difference is only that the angles of such edges are computed after the transformation M _1 
is applied to them. 
Thus, we reduced 

BondVol(T, r) = BondVol(M, MT, Mr) = BondVol(M, MT, e 3 ). 



3.1.3 Reduction to Truncated Prism 

Denote the vertices of T by A, B, C, and D; and choose their ordering to have a positive orientation 
in 3D: i.e., so that the vectors AB, AC, and AD form a positively orientated basis. 

The tetrahedron T can be represented as an oriented sum of truncated prisms, which can be 
rigorously expressed with characteristic functions: 

X T = -o(B'C'D') Xp{BCD) + o(A'C'D') Xp{ACD) - o(A'B'D') Xp(ABD) + o(A> B'C) Xp{ABC) , (3.8) 

where by we denote a projection of a point or a vector on the xy-plane (i.e., on the plane 
orthogonal to 63), P(XYZ) is a truncated prism with vertices X, Y and Z and their projections 
X', Y', Z' (see Fig. 



3.1 



for an illustration), and o(XYZ) G {—1,0,1} is an orientation of three 
points X, Y, Z G defined to be zero if X, Y, Z lie on the same line and to be equal to the 



XZ 



otherwise. 



orientation of the basis XY. 

The lower-dimensional version of (3.8) (i.e., splitting of a triangle into trapezia) is illustrated 
in Fig. 2(a) and the proof for an arbitrary dimension is given in Appendix |A| 

For convenience, we assume that all four points, A, B, C, and D, lie above the xy-plane 
(otherwise some prisms may be ill-defined). One, obviously, can always shift T upwards to satisfy 
this requirement. This, however, is not required with an appropriate generalization of Xl 
when T is not entirely above the xy-plane; refer to Appendix [A] for more details. 
Thus, we reduced computing BondVol(T, r) for a tetrahedron T to computing 



"-P(XYZ) 



BondVol(M,T, r) 



- o(B'C' D')BondVol(M,P (BCD),r) 
+ o(A'C'D')BondVol(M,P{ACD),r) 

- o(A / S / D / ) Bon dVol(M,P(ylSD),r) 
+ o(A , J B , C / )BondVol(M,P(,4 J BC7),r). 
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Figure 3.1: A truncated prism P(XYZ) formed by three vertices X , Y, Z and their projections, 
X' , Y' and Z', on the xy-plane. 



3.1.4 Reduction to Sums over Triangles 

We have M G Z 3x3 , det M = 1, r = e^, P = P(ABC) is a truncated prism, the points A, B, and 
C are above the xy-plane; and we need to compute BondVol(M, P, r). In what follows we will use 
the notation A(XYZ) for a triangle with three vertices X, Y, Z G M 2 . 

We can assume that the plane ABC is not parallel to the z axis: otherwise P(ABC) is degenerate 
and BondVol(M, P, r) = (since then Xp = 0). Hence, let z = c\x + C2y + C3 be an equation of the 
plane ABC. 



We will split all the bonds, as we did in (3.5), into the classes + r), £ G (cf. (3.2)), each 
class defined by x = ie\ + je% , and i,j G Z. That is, we express 



BondVol(P, r) 



EE/ ^db 



/OO r 
x¥(i,j,z)dz= \ 



cii+C2j+c:j 



Xp(i,j,z)dz. 



For X3 between and c\i + c%j + C3, Xp(*> Js z ) is equal to 1, |, or, respectively, depending on 
whether is in the interior of A = A(A' B'C), on its edge, or, respectively, on its vertex. The 
value a is determined by the edges, v and w, sharing the respective vertex: a is equal to the angle 
between the plane formed by M -1 n and M _1 e3 and the plane formed by M~ 1 w and M _1 e3. The 
latter is equal to the angle between M -1 ?; x M _1 e3 and M -1 ?/; x M _1 e3. 
Thus, if we define, for a polygon Scl 2 , its characteristic function 



1 if G interior of S 
\ if G edge of S 

£ if G vertex of S sharing edges v, w, where 

a is the angle between M~ 1 v x M _1 e3 and M w x M -1 e3, 



then 



Cli+C2j+C3 



Xp(h3,z)dz = (cii + c 2 j + C3)x™( A , B , c/) (i,j), 



~M 
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Figure 3.2: On the left: A triangle A(ABC) and the projections of its vertices on the x-axis. 
The triangle can thus be represented as an oriented sum of three trapezia, Trap(AB), Trap(-BC), 
Trap(CA). One can notice that the area under A(ABC) is counted once with minus (for Trap(AB)) 
and then with plus (for Trap(i?C) and Trap(CM)). On the right: splitting of a trapezoid Trap(AB) 
into a right triangle A{ABD) and a rectangle Trap(AD). 



and hence 



o(A'B'C')BondVol(P,r) = o(A'B'C') ]T {c x i + c 2 j + c 3 )x™ {A , B , cl) (i,j) 



'■.it 



-■ S' A , B/cl (M;c 1 ,c 2 ,c 3 ). (3.9) 
We thus reduced the problem to computing the sum S' A , B , C , (M; 01,02,03) over a triangle A(A'B'C'). 



3.1.5 Reduction to a Sum over Right Triangles 

Let A = A(ABC) and let A', B' , and C be the projections of A, B, C on x-axis (i.e., on e\), as 
illustrated on Fig. 2(a)| 
Then 



o(ABC)x^ 



(ABC) 



-o(BC)x™ MBC) + o(AC)x^ MAC) - o(AB)xt M AB), 



(3.10) 



where Tvap(XY) is a trapezoid with vertices X, Y, X' , and Y' , by •' we denote a projection on the 
x-axis (i.e., on ei), and o(XY) is the orientation of the two points X and Y on the x-axis defined 
as sgn(Ai • ei). 

The formula (3.10) is quite intuitive: Indeed, as seen on Fig. 2(a)[ the area under the triangle 
will be counted with the minus sign when evaluating — °(^4^)Xxrap(AB) an< ^ w ^ n the plus sign when 
evaluating -o(-BC)XTra P (BC) + (^^)^Trap(AC*)- The P r001 * °f ( |3.10| ) is given in Appendix 



A trapezoid Trap(Al?) can further be split into a right triangle and a rectangle (see an liTustra- 



tion on Fig. 2(b) ): 



o(AB)xt M AB) = o{ABD)xt {A BD) + o(AD)x% { ad), 
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Figure 3.3: A right triangle A (ABC) (left) and its shifted copy A (ABC) together with its rotated 
copy A(Ct>i) (right). 

where D := (Bi,A 2 ). 

Thus, we reduced our problem to two problems: (1) computing S' ABC (M; c±, c 2 , C3) where AB 
and BC are parallel to x and y axes respectively, and (2) computing 

°{ AD )X^p{AD)(hj)( C ^ + C2J ' + C3 )' 

where AD is parallel to the x-axis. The latter can be computed analytically using the fact that the 
function °( A D)x^j. ap ^ AD ^{i, j) ^ s sy mme tric w.r.t. rotation by the arc length ir around the center O 
of the rectangle Trap(AD): 

o(AD)x^ ap{AD) {iJ){ Cl i + c 2 j + c 3 ) = (D 1 - A l )A 2 {c x i + c 2 j + c 3 )| (iJ)=c ,. (3.11) 
3.2 Computing the Sum over a Right Triangle 

It remains to develop an algorithm of computing S' ABC (M; 01,02,03) for AB and BC parallel to x 
and y axes respectively, where S' is defined by (|3.9l). 



Let A = (Ai,A 2 ), B = A + be±, C = B + ae 2 (see Fig. 3(a) ). We assume that both a and b are 
nonzero (otherwise A (ABC) is degenerate and S' ABC (M; c%, c 2 , C3) is zero). 

We shift the points A, B, and C so that A coincides with the origin (see Fig. |3(b)[ ). That is, 
we introduce the points A = (0,0), B = (6,0), C = (b,a), and change the variables of summation 
i -> i - Ai, j -> j - A 2 : 

S'ABc( M ; c hC2,C 3 ) = o(ABC) ^ XA(ABO)( i + A ^j + A 2)(cii + C 2 j + C4) 

=■ 5 , ^^(M;ci,c 2 ,c 4 ), 

where 04 = 03 + c\A\ + C2A- 

Second, notice that since A(ABC) and its copy rotated by it, A(CDA) (see an illustration on 
Fig. |3(bj| ), together compose a rectangle, we can use (3.11) and express 



5^(M;0,0,C4) = i(^ M (M;0,0, C 4) + ^ M (M;0,0,c 4 )) = \abc A . 
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Thus, it remains to compute S^^(M; C\, c 2 , 0). 

Third, notice that we can reduce it to the case a, b > by doing reflections w.r.t. the axes (e.g., 
reflection around ei-axis corresponds to changing a — > —a, c\ — > —c±, o(ABC) — > —o(ABC)). 
Hence we assume that a, b > and therefore o(ABC) = o(ABC) = 1. 

Last, notice that the function Xa(ABC)^'-^ can ^ e described as follows: 



A 



A(ABC) 



0<i<b, < j < § i 
0<i<b, j = 
< i < b, j = ji 



i = b, < j < a 
i = 0, j = 
i = b, j = 
i = b, j = a 
otherwise, 



2tt 
2tt 



with 



a 

P 
7 



ang(M _1 (6ei) x M _1 e 3 , M _1 (6ei + ae 2 ) x M _1 e 3 ) 
ang(M~ 1 (6ei) x M _1 e 3 , M _1 (ae 2 ) x M _1 e 3 ) 
ang(M _1 (ae 2 ) x M _1 e 3 , M~ l (bei + ae 2 ) x M _1 e 3 ), 



where ang(v, w) := arccos(t> • w ) for u,v £ M 3 . 
Thus, 

6-1 L?J 

S' A6d (M;a,c 2 ,0) = ^^(cii + c 2 j) 

i=l j=l 



6-1 



\ M + C2f) + ^(ci6 + c 2 j) + ^( Cl z + c 2 0) 

3=1 



0<i<b 
j£gcd(a,6)Z 

+ E ~~' M 

(*j)e{i4,B,<3} 



i=l 



Each of the terms in the second line is a sum of an arithmetic progression and can be expressed 
analytically. 

Thus, it remains to compute 

6-1 Lf J 6-1 Lf J 

1=1 3=1 i=0 j=l 

In what follows we will use the following two standard identities that can be easily proved by 
induction: 



n-l 
i=0 



n(n — 1) 



n— 1 



£° 2 = 

i=0 



2 n(n-l)(2n-l) 



(3.12) 
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Using the second identity we can transform 
6-1 L^J 6-1 . 6-1 



i=0 j=l 



i=0 



i=Q 



ai 1 



E E < = E4tJ = £«'lf - 6 (M mod 6) J = 6 (6 ^ : ),,i2/ ' - 1 1 - 



where 



6-1 



i=0 



^,6 := E ( ai mod h ^ 



(3.13) 



and 



6-1 Vf\ 

EE^ 

t=0 j=l 



6-1 



1 | 

6~ 



a/ 



+ i 



- >^ I — — T\ ai mod on II + — — T\ ai m od oj 

Zi \ (J (J / \ U (J 

i=Q 
,6-1. 

1 v-^ az / cm 

2Et( 1 + it 

i=0 



6-1 



1 6-1 1 6-1 ai I — 

— ~^2(ai mod b) — - "T~( a ^ m od b) + y~](ai mod b) z 

.11 (J (J /.it 

i=0 



2b 



i=0 



i=0 



The first, second, and fourth sum can be computed analytically: 

6-1 



Eai ( at 
i=0 



1)(36- a + 2ab) 
66 : 



6-1 gcd(a,b) 

(ai mod b) = gcd(a, b) ^ gcd(a, b) i = 



i=0 
6-1 



i=0 



» _i 

gcd(a,6) 



y^(ai mod b) 2 = gcd(a, b) (gcd(a, b) ■ 



i=Q 



i=0 



2 _ b(6-gcd(a,fe))(26-gcd(a,b)) 
J 6 



(3.14) 



(3.15) 



(3.16) 



and the third sum again reduces to (3.13): 



1 6-1 ' 

lEf( ai mod 6 ) = I s - 



b^ b 

i=0 



b- 



Here the identity (p. 14) follows from the standard sums (3.12); and to compute (3.15) and (3.16) 



we used (i) the fact that the numbers (ai mod b), i = 0, 1, . . . , b — 1 are essentially the numbers 
{0, gcd(a, b), . . . , ( g C d(a 6) — -0 § ca -( a ' k)} repeated gcd(a, b) times each, and (ii) again the standard 
sums (3.12). 



3.3 Computing The Sum S a , 



For computing the sum (3.13) we propose a Euclidean-like algorithm with complexity O (log (a +6)), 



which consists in iterative reducing the (a, 6)-problem to (b mod a, a) problem. 
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Using the following identity [5, p. 85]: 



a-l 



[ ax \ = 

j=0 



X + 



Vx 6 



with x = i/b (we assume b ^ 0, since 6 = is a trivial case), express 

a-l 



-(az mod b) = — 



Hence transform 
6-1 . 



i=0 



6-1 



S a ,b = ^2l( ai mod b ) = ^2' 



i=0 



V 



a-l 

E 

3=0 



a i 
J 



% J 
b a 



a i 
J 



E 

3=0 



6-1 



1 J 

b a 



6-1 a-l 



E» 2 -EE> 

i=0 3=0 



i=0 



o a 



=: Si — S 2 . 



The first sum is trivial: Si = g(6 — l)a(26 — 1). 

For computing the second sum, notice that \i + equals to 1 if £ + ^ > 1 and otherwise, 

and 

bj 



\ + 3 ->l ^ i>b-— ^ i>b 
b a a 



Hence 
S 2 



6-1 a-l 

EE- 

i=0 i=0 



1 J 

b a 



a-l 6-1 



a-l 



E E > = E 

3=0 i=b _ ^Mj 3=0 





(2b - 1 - 




) 


a 




a 





- ( — (bj mod a) J [2b — 1 H — (6j mod a N 

3=0 v 7 v 



a-l 



1 ^ (b{2b- 1) 

2 ^ 

3=0 



3 n3 

a 1 



b 2 . 2 26-1 



(bj mod a) H — ^ ,7(6,7 mod a) ^-(6,7 mod a) 



We next compute individual sums in S2. 

Using (3.15) and (3.16), evaluate the individual sums in S 2 : 

(a-l)(2a-l)6 2 



S 2 = ~(a- 1)6(26-1) 



12a 



t (I -2b) (a-gcd(a,6))a | 6 (a - gcd(a, 6)) (2a - gcd(a, 6)) 

+ 2a 2 a 6 ' a 12a 

Substituting this back to S a ,b and collecting the terms yields 

3a 2 6 + 3a6 2 + a 2 - 3a6 + 6 2 - 6a6 gcd(a, 6) + gcd(a, 6) 2 6 

'-'a, 6 — 77i — &b,a- 

12a a 

Finally, notice that Sb >a = Sf, mo da,a- We thus managed to reduce the problem from (a, 6) to 
(6 mod a, a), which yields the C(log(a + 6)) algorithm. 
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3.4 Complexity of the Algorithm 



Proposition 2. The number of operations of the algorithm described in Sections 3.1-3.3 is at most 
e>(log(diam(T)) + log(|r|)). 



Proof. Evidently, \r\ is not increased as a result of reducing to gcd(ri, r 2 , r%) = 1 (Section |3.1,1 ), 



In reducing to r = e% (Section 3.1.2), we can choose < c\ 2 < | r*2 1 , < du < \fi\, < C3 < | r*3 1 



< d% < gcd(|ri|, |j*2|). Then, up to a constant factor, ||M|| can be estimated as 

M+H 

gcd(|n|,|r 2 |) 



M||oo < max(|r 3 | + gcd(|n|, |r 2 |), 1, c 3 + d 3 ) max (c 12 + d 12 , 'gnMrll) ■ 1 ) ~ l r | 2 - 



And hence diam(MT) < r 2 diam(T). 



The triangles considered in Section 3.1.4 are the bases of the truncated prisms of Section 3.1.3 



which are projections of faces of MT and hence the diameter of each A in Section 3.1.4 is at 
most 0(r 2 diam(T)). Similarly one can deduce that both a,b in each sum S a b are of the or- 



der £?(r 2 diam(T)) and hence the Euclid-like algorithm of Section 3.3 takes C( log(r 2 diam(T))) 
0[ log(diam(T)) + log |r|) operations. □ 

4 Numerical Tests 

Numerical tests were conducted in order to numerically study the accuracy and stability of the 
method. 

The atoms interact with Lennard-Jones potential 4>{z) = — 2|z| -6 + |z| -12 under which the FCC 
lattice is stable. The cut-off radius is chosen to be 3.2. 

The lattice vectors are chosen as a\ = (0>7y3>7^|)> a 2 = (^/f'^'^/j)' a 3 = (^75; 7^3 '^)- The 

reference lattice C \ £d consists of a crystal whose atoms formed a cube with the side 2^/2N 

(N = 8, 16). The Dirichlet-type boundary conditions are imposed by introducing additional atoms 

with fixed positions, £d- A single atom at the origin is removed thus forming a vacancy defect. In 

total, #(£ \ £d) = 125022 (unconstrained) atoms are in the atomistic system. 

/1 0.01 0.02 \ 

A macroscopic uniform deformation with F = I 1 0.015 I is applied to the constrained 

Vo 1 / 

atoms £d- 

The atomistic region Q a is chosen as a smaller cube with the side 2\/2K, K = 2, 3, ... ,11 (see 



an illustration on Fig. 4.1 ). A quasiradial mesh with the mesh size h = 1 near the A/C interface is 



chosen in accordance with the optimal choice of meshes in 2D [13J. 
4.1 Accuracy Test 

In the numerical tests, the exact and the approximate solutions were computed using Newton's 
method of solving the equilibrium equations, with the initial guess being an undeformed configu- 
ration. 



The results of the computations are shown in Fig. 4.2 The difference in the W 1,00 -norm between 
the approximate and the exact deformation is plotted on the left and the difference between the 
energies \E h {u h ) — E{u)\ is plotted on the right. One can observe at least a third order convergence 
in the W 1,00 -norm and close to fifth order convergence for the energy. 
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Figure 4.1: Geometry of the a/c coupling: a triangulated continuum region and atoms in the 
atomistic region, for K = 2. 




Figure 4.2: Results of computations. The W 1,00 -error (left) and error in energy (right) are plotted 
against the refinement level K = 2, 3, . . . , 11. One can see a third order convergence in the W 1,0 °- 
norm and close to fifth order convergence for the energy. 
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Figure 4.3: Stability regions for a Bravais lattice. The solid line corresponds to the exact stability 
region of the infinite lattice, the dashed line corresponds to the coupling with N = 16 and K = 8, 
and the dotted line corresponds to the coupling with N = 8 and K = 4. The results suggest that 
the coupled A/C system does not lose stability earlier than the original atomistic system. 



4.2 Stability Tests for a Bravais Lattice 

We also conduced stability tests for a Bravais lattice (i.e., with no defects) to verify that the stability 
region of the proposed coupling E h is not smaller than the stability region of the atomistic model 
E. This is crucial in numerically studying defects: it must be made sure that the onset of instability 
occurs due to motion of a defect but not due to artifacts of the coupling. 

We take a Bravais lattice C\Cd similar to the previously described but with no removed atoms. 

/l + t 0.05 0.02\ 

The macroscopic uniform deformation gradient F = I o l + s o.oi I is applied to the constrained 

V o o l / 

atoms Cd. The two parameters, t and s, were varied in the range between —0.2 and 0.2. 

We computed stability regions for the coupling E h and compared it to the stability region of the 
exact atomistic model on an infinite lattice. A stability region is defined as the set of parameters 
(t, s) for which the model is stable. Stability of E h was determined by numerically testing whether 
the Hessian 5 2 E h is positive definite. The stability of the atomistic model was computed analytically 
with the help of Fourier transform [6] . 



The stability regions are plotted in Fig. 4.3 The solid line corresponds to the exact stability 
region of the infinite lattice, the dashed line corresponds to the coupling with N = 16 and K = 8, 
and the dotted line corresponds to the coupling with N = 8 and K = 4. One can see that stability 
regions of the coupling strictly contain the exact stability region and seem to approach it as N and 
K increase, which is the desired behavior of the A/C coupling. 
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5 Discussion and Conclusion 

In the present paper the consistent A/C coupling [16| has been extended to 3D for two-body 
potentials. The proposed method couples the atomistic equations with the modified Cauchy-Born 
continuum model. The continuum energy of the modified model can be evaluated efficiently, as 
discussed in Section [3j Although the stability of such a modified continuum model has not been 
studied analytically in the existing literature, the numerical tests suggest that the proposed coupling 
is stable whenever the atomistic model is stable. The numerical tests also confirm convergence of 
the proposed coupling to the exact solution. 

The major challenge yet to be solved is an extension of the present method to many-body 
interaction. A one-dimensional consistent coupling for such an interaction exists [9], however it 
does not seem obvious how to define continuum approximations to many-body bond energies in 
many dimensions. 
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A Splitting a Simplex into Truncated Prisms 
A.l Auxiliary Definitions 

We define an orientation of points X^\ . . . , X^ d+1 ^ € R rf as an orientation of the basis (X^ 
X<n,...,xW)-XWy,Le. 



/ Y^ Y^ Y^ 

A l — A l A 2 — A 2 



o d (I (1) ,..,X (W) ) :=sgnodet 



xf - xP 



y(3) y(l) 
A 2 ~~ A 2 



\ A 1 _A 1 A 2 _A 2 



xf-xf 



4 d+l, -4 1 V 



, (A.l) 



o d (xm,...,X^)€ {-1,0,1}. 

Next, we define the characteristic function of a truncated prism P = P(X^,...,X^) for d 
points X^\ . . . , X^ £ M d . If the plane through these d points is perpendicular to the hyperplane 
Xd = then we let x P '■= 0- Otherwise let x d = Yli=i a i x i be an equation of such a plane and 
define x P G Li( Rd )> 



X 



p '" 



'1 < x d < Eti 1 and £ e conv((X( 1 ))', . . . , (X^)'), 
-1 > x d > E-=i CHXi and £ G conv((A'( 1 )) / , . . . , {X^)% 
otherwise, 



and 



Xp {x) := lim 1 / Xp(O d C, 
^ p^o \Bp[x)\ J Bp ( x ) F 



(A.2) 



(A.3) 
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where by we denote an orthogonal projection on the hyperplane x d = 0. 

Notice that X P = X P almost everywhere (x p is defined pointwise in M rf and may take interme- 



diate values conforming with the definition (2.6)). Also notice that x P is a characteristic function 
of a polytope P in the sense of (2.6) only when the face conv(X^ 1 \ . . . , XW) is "high enough" (i.e., 



has a large enough a^-coordinate) . Otherwise, writing P = P(X^\ . . . ,X^) is formal and does 
not refer to a proper polytope in R . 



A. 2 Formulation of the Result 

Let A^\ A^ d+l "> £ R d so that o d (A^ , . . . , A^ d+1 '>) = 1. Denote 



the simplex T 
the truncated prisms P^ 
their orientations 



and the faces of T, 



i?( fc ) 



conv 

od-iaAtty, . . . , (A^y, (A( k +vy, . .'. , (A^y), 



conv 



(aw,...,^- 1 ),^ 1 ),...,^ 1 )), 



where k = 1, . . . , d+ 1. Here we identify points on the hyperplane x<2 = with points in M^ -1 when 
computing orientations via (A.l) 
We prove 



Proposition 3. 



d f(-^ k) x pw 

k=l 



(A.4) 



pointwise in E , where x T is defined by (2.6). 



A.3 Proof 



In the following two lemmas we prove that (A.4) holds almost everywhere. 



i a (k) 

Lemma 1. The relation (A.4) holds almost everywhere in M. if A , > for all k = 1, . . . , d + 1 



(i.e., ifT is entirely above the hyperplane x d = 0). 



Proof. Since all vertices A^ lie above the hyperplane x d = 0, the function x P defined by (A. 2) 



and (A.3) equals 1 or almost everywhere and hence is the characteristic function of a proper 



truncated prism P. 

Fix an arbitrary test function / £ C°°(]R c( ), and let 



g(xi, . . . , xa-uXd) := e d / f(x 1 , x d -i,£)d£ 

Jo 



so that / = divg (recall that e±, . . . , e d is the canonical basis of M. d ). Multiply (A.4) by / and apply 
the divergence theorem: 



d+l 



g ■ n T d7 



OT 



fc=i 



/ G (fc) 
dp( k ) 



w 5 -n p(fe) d 7 , 



(A.5) 



where n* is the outward normal vector (to T or P^ respectively). 
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Next, notice that g = on the base of the truncated prims (i.e., for x d = 0) by the definition of 
g and g • n p< ~ k) = on the sides of the truncated prism (i.e., below F^). Hence, in both parts of 
the relation (A. 5) we have the sum over faces F^ k > of the integrals of g ■ n* and it only remains to 
verify that n T = — {—l) k o^n p k for all k. 

To prove that, fix k S {1, . . . , d + 1} such that is not degenerate (otherwise the statement 
is trivial since n p ^ k = and n T = on that face of T). Notice that n p< ~ k) > (i.e., the vector n P (k) 
points upwards). The agreement of orientation of n T and o^ k ' follows from the following chain of 
statements: 



(-1)* 



& o d {A {1) A (k ' l) , A^ , A® + ed ) = (_ if 

o d {A (1) , + e d , A^ A^) = 1 

44> e d is an inward vector w.r.t. 
O n d < 0, 

where £ is any integer between 1 and d + 1 different from k. Here in the first step we used 
the expansion of determinant by minors and in the third step the following fact: Since the basis 
(A^ - AW, . . . , - AW) is positively oriented, a vector v is an inward (or, resp., outward) 

vector w.r.t. if and only if the orientation of the basis (A^-A^, A^ k ~^ - A&\ v, A^ k+l ^ - 
A^\ A^ d+l ^ - AW) is positive (or, resp., negative). □ 



Lemma 2. The relation (A.4) holds almost everywhere in 



Proof. In view of the previous lemma, we only need to show that both sides of (A.4) are invariant 
w.r.t. Sd, a dilatation in x d by distance D. More precisely, we need to notice that SdX^ = X c <r 
almost everywhere (which follows directly from the definition of x) an d prove that 



d+l d+l 

Y,(-^ k o {k) s D o Xp(k) 

k=l k=l 



(fc) 



x 



(A.6) 



Proof of (A.6) is based on noticing that 

X S D P(XW,...,X( d )) = Sd ° X P(XW,...,X( d )) + X P(De d +(XWy ,...,{De d +xwyy 

i.e., that x °f a shifted truncated prims equals to shifted x °f a truncated prism plus x °f a prism 
of height D with the base formed by prjections (X^)' , . . . , (X")'. This can be verified by fixing 
x £ conv((X^ 1 -') / , . . . , (X^)') (for x elsewhere the statement is trivial) and considering three cases: 
Yli=x a i x i is (i) less than min(0,D), (ii) between min(0,.D) and max(0,D), and (iii) greater than 
max(0, D). 



We now take the difference between the left-hand side and the right-hand side of (A.6): 



d+l 



£(-i: 



k n (k) 



Sd°X 



P (k) 



k=i 



£(-i) fc (fe) 

k=l 



x s D p( k '> 



d+l 



x s D (F( k )y 



(A.7) 
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Figure B.l: An alternative splitting of A (ABC) into three right triangles and a rectangle. 



and notice that this difference depends only on projections of faces, (F^)' , but not on d-th coordi- 
nate of the points AW, . . . , A^ d+1 \ Hence we can again shift these points so that Lemma [l] applies 



to both T and SpT and hence conclude that the difference (A.7) equals to Sp ° X 
almost everywhere. 

We now finalize the proof of Proposition [3| 



X 



S D T 



□ 



Proof of Proposition^ Denote f(x) = x T i x ) + Ylk=i(~ l) fe °^ X P (k)( x )- From (2.6) and (A. 3) we 
have 

/(x) = lim— L- / /(£)d£. (A.8) 



p->o \B p (x)\ J Bp{x) 



Due to lemma[2], /(£) = for almost all £, hence the right-hand side of (A.8) is zero, hence f{x) = 
for all x eR d . □ 



B A Matlab implementation of computation of BondVol(T, r) 

The code given in this appendix closely follows the algorithm outlined in Section [3] except for the 
following optimization. 



Instead of the splitting (3.10) we introduce three additional points D = (B\, A2), E = (Ci, A2) 



F = (CijBz) (see Fig. B.l) and represent 

1 

\(ABC) 



o(ABC)x™ 



- o{ADB)x^ (ADB) 
+ o(AEC)xYL [AEC) 
+ o{FBC)x*t (FBC) 
+ o{EDB)x™ {EDBF) , 

where U\(EDBF) is the rectangle EDBF. The contribution of the latter is computed by a formula 



analogous to (3.11). 
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The Matlab Code 

function BondVol = BondVol_tetrahedron(v, r) 

% BondVol (v,r) of the tetrahedron with vertices v 

7, change the coordinates so that r = (0,0,1) 
[gcdl2, cl2, dl2] = gcd(r(l) ,r(2)) ; 
[gcd3, c3, d3] = gcd(gcdl2,r(3)) ; 
if (gcd3~=l) 

7, reduce to the case gcd3=l 

r = r/gcd3; 

gcdl2 = gcdl2 /gcd3; 

end 

M = eye (3) ; 

if (r(l) ~=0 I I r(2) ~=0 I I r(3) == 0) 

M= [cl2,dl2,0;r(2)/gcdl2, -r(l)/gcdl2, 0;0,0,1]; 
M = [-r(3) ,0,gcdl2;0,l,0;c3,0,d3]*M; 

r = M*r; 
v = M*v; 

end 

invM = inv(M) ; 

BondVol = BondVol_prism(v(: , [1;2;3]) , invM) + ... 
BondVol_prism(v(: , [4;3;2]), invM) + ... 
BondVol_prism(v(: , [4;2;1]), invM) + ... 
BondVol_prism(v( : , [1;3;4] ) , invM) ; 

end 

function BondVol = BondVol_prism(v, invM) 
7. BondVol (v,r) of the truncated prism, 

7o formed by three vertices v and their projections on the XY plane; 
'/, r is assumed to be (0,0,1); 

7 invM*r and invM*v are the original positions. 

7. shift the triangle (and the prism) so that vl = (0,0,*) 
v([l 2],2)=v([l 2],2)-v([l 2],1); 
v([l 2],3)=v([l 2],3)-v([l 2],1); 
v([l 2] ,1) = [0;0] ; 

if (round (det([[l;l;l] , v([l 2] ,:)']) )==0) 
7o if the prism is degenerate 
BondVol = 0; 

else 

7. find coefficients [c4; cl; c2] of the function cl*i + c2*j + c4 that we sum 
c4clc2 = [[1;1;1], v([l 2] , : ) ' ] \v(3 , : ) ' ; 

7. reduce to integration over a triangle formed by (0,0), (v2x, v2y) , (v3x, v3y) , 
7. which is further reduced to integration over three right triangles and 
7. one rectrangle 

BondVol = + right_triangle_sum(c4clc2, v([l 2], 2), invM) ... 
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- right_triangle_sum(c4clc2, v([l 2], 3), invM) ... 

- right_triangle_sum(c4clc2 + ... 

[1;0;0]*([0 v([l 2] ,3) '] *c4clc2) , v([l 2],2)-v([l 2], 3), invM) ... 
+ prod(v([l 2] ,3)-[v(l,2);0])* ([1 0.5*(v([l 2] ,3) + [v(l , 2) ; 0] ) >] * c4clc2) ; 

end 

end 

function ans = right_triangle_sum(c4clc2 , pt, invM) 

°/ sum c4 + cl x + c2 y over a right triangle (0,0), (pt(l),0), pt 

b = pt(l) ; 7, x-side 
a = pt(2) ; 7, y- side 

if (a==0 I I b==0) 

ans = 0; return; 

end 

7o reduce to a>0 and b>0 

orientation = sign(a) *sign(b) ; 

c4clc2 = c4clc2 .* [1 ; sign(b) ; sign(a)] ; 

invM = invM*diag( [sign(b) , sign(a), 1]); 

a=abs(a); b=abs(b); 

7, sum the constant term 
ans = l/2*a*b*c4clc2(l) ; 
clc2 = c4clc2(2:3) ; 

7o sum the linear terms, using reduction to Sab 
gcdab = gcd(a,b) ; 
Sab_ans=Sab(a,b,gcdab) ; 

ans = ans + c4clc2(2) * (l/6*a* (b-1) * (2*b-l) - Sab_ans) ; 
ans = ans + c4clc2(3) * (1/4* (a-1) * (b-1) + 1/4* (gcdab-1) ... 

+ l/12*(a~2)/b*(b-l)*(2*b-l) + 1/12/b* (b-gcdab) * (2*b-gcdab) - a/b*Sab_ans) ; 

7o sum over sides: 

ans = ans + 0.5*(a-l)* [b a/2]*clc2; 
ans = ans + 0.5*(b-l)* [b/2 0]*clc2; 

7o sum over the hypotenuse (we subtract half the contribution) 
ans = ans - . 5* (gcdab-1) * [b/2 a/2]*clc2; 

7o sum over vertices 

vl = cross(invM(: ,1) ,invM(: ,3)) ; 7, invM(:,l) == invM*[l;0;0] 

v2 = cross(invM(: ,2) ,invM(: ,3)) ; 

v3 = -cross(invM*[b;a;0] ,invM(: ,3)) ; 

i = 0; j = 0; 

ans = ans + acos(dot(vl,-v3)/norm(vl)/norm(v3))/(2*pi) * [i j]*clc2; 
i = b; j = a; 

ans = ans + acos(dot(-v2,v3)/norm(v2)/norm(v3))/(2*pi) * [i j]*clc2; 
i = b; j = 0; 
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ans = ans + acos(dot(-vl,v2)/norm(vl)/norm(v2))/(2*pi) * [i j]*clc2; 
ans = real(ans) * orientation; 

end 

function ans = Sab(a,b, gcdab) 

% sums i/b (a*i mod b) over i=0..b-l 

if (b==0) 

ans=0; return; 

end 

result = 0; multiplier = 1; 

while(a ~= 0) 

result = result + multiplier * ... 

(3*b*a~2 + 3*b~2*a + a~2 - 3*a*b + b~2 - 6*a*b*gcdab + gcdab~2)/(12*a) ; 
multiplier = multiplier * (-b/a) ; 
old_b=b; b=a; a=mod(old_b,a) ; 

end 

end 
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